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1. Introduction 

CM 

In this work we discuss a numerical method for studying singularities (blow-ups) in the Keller- 
Segel model (K-S) of chemotaxis 1 28 1 and in the systems with similar mathematical structure, such 
as the McKean-Vlasov (McK-V) model of charged interacting particles EDI . (Patlack derived simi- 
lar equations nearly two decades prior to Keller and Segel [29 j.) We present illustrative numerical 
i '~ l i examples and provide analytical insights into formation and interaction of singularities from the 

point of view of the associated stochastic processes. We concentrate on the details of numerical 
modeling and several specific aspects of mathematical analysis. An interested reader should con- 
sult the exhaustive reviews for discussions of physical and biochemical phenomena involved in 
chemotaxis 



> 
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Numerical treatment of equations with singular solutions is challenging because it requires a 
reliable and efficient approximation of singular functions or even distributions. For example, the 
conventional finite-difference or spectral methods are not fit for dealing with ^-functions arising in 
the K-S model. While most finite-element and discontinuous Galerkin methods are efficient and 
accurate as long as the solutions remain sufficiently smooth, in the vicinity of singularities they 
• • struggle with positivity issues, unphysical oscillations, and related problems ll20llT9llT7llT8ll34ll33l . 

(A conservative upwind finite-element method introduced by Saito for the elliptic K-S model over- 
comes some of these issues |30J.) Generally, unless some regularizations are introduced, these 
methods cannot propagate solutions past the moment of blow-up altogether. The method pre- 
sented in this work avoids such difficulties because the singular field is modeled by an ensemble of 
interacting particles. The singularities manifest themselves as particle aggregates, harmless from 
the numerical standpoint. Essentially, this approach to chemotaxis reverses the continuous PDE de- 
scription of ensembles of interacting particles. Haskovec and Schmeiser have recently introduced a 
method utilizing similar ideas I12T1 I22I. Their method, however, is exclusively particle-based and is 
limited to one specific variety of the elliptic K-S model in the entire plane. Our method employs the 
grid and particle representations of the fields simultaneously and is suitable for both elliptic and 
parabolic K-S equations in arbitrary domains. It is also applicable to a wider range of kinetic prob- 
lems. (Other technical differences are discussed where appropriate.) Conceptually, our method 
is more in the spirit of the Particle In Cell (PIC) methods widely used, e.g., in plasma physics 
ll2^IT4l l5ll39l[Tl l23ll4TII . Another interesting numerical method has been recently implemented for 
the one-dimensional McKean-Vlasov system |6|. It uses that the McK-V equation describes steepest 
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descent dynamics in the Wasserstein-2 space of probability measures and is suitable for studying 
singular (measure-valued) solutions. 

The usefulness of the particle approach extends beyond numerical methods. Analysis of the 
stochastic processes underlying the K-S PDEs yields insights into formation and interaction of sin- 
gularities. In particular, we derive an expression for the critical mass required to create or sustain 
a singularity and relate the non-uniqueness of solutions of the K-S equations to the underlying dif- 
fusion process. Even though our treatment is rather informal, it sheds new light on the mechanism 
of blow-ups. Some of the results presented in this paper have been proven using different methods 
and under various assumptions. Much effort has been directed towards obtaining the exact value 
of the critical mass |15, 8, 9| [TTJ [4l . Velazquez introduced measure-valued solutions and extended 
the K-S equations beyond the point of blow-up (33120/ see a ls° ED- These works together with 
the study of the K-S equations as a hydrodynamic limit of interacting particle systems by Stevens 
1 32] link the particle and PDE descriptions of chemotaxis. We also mention a rigorous analysis of 
phase transitions in the McK-V system [13] (it only covers sufficiently regular interaction kernels 
not relevant to the K-S model) and an investigation of singularity formation in several related ag- 
gregation models |2. 3, 10] (these models do not have diffusive terms and correspond to one special 
case of the K-S model). 

The Keller-Segel model is prescribed by the following system of PDEs: 

d t p(x,t) = V-{pVp - XpVc), (la) 

ol d t c(x,t) = Ac - k 2 c + p. (lb) 

The function p(x) is the density of active particles (bacteria), c(x) is the concentration of chemoat- 
tractant. For numerical simulations we use Neumann (no flux) boundary conditions, 

d n c(x) = d„p(x) = for xedD, (2) 

in a two-dimensional square domain O = (0, L) 2 . We assume that p(x) integrates to M over the 
entire domain; M is the total mass of the bacteria. By rescaling the equations, we can always achieve 
that a = or a = 1. In the former case the model is called elliptic, in the latter — parabolic. The 
parameters p, Xi and k are constants. (The numerical method, however, may be readily extended 
to a more general class of equations, see Section|2]for details.) 

Physically, these equations describe the following phenomenon. The chemoattractant spreads 
diffusively and decays with rate k 2 ; it is also produced by the bacteria with rate 1. (In the ellip- 
tic case these rates are infinite, i.e., the chemoattractant "thermalizes" infinitely fast.) The bacte- 
ria diffuse with mobility p and also drift in the direction of the gradient of concentration of the 
chemoattractant with velocity x\ Vc|; x 1S called chemosensitivity. 

In this work we are interested in singular solutions to the K-S equations. A typical situation 
is illustrated in Figure [I] where a snapshot of the concentration field c(x) is displayed. The peaks 
correspond to ^-function-type singularities of the particle density p (x) . (Displaying bacteria density 
itself is not very illustrative due to its highly singular nature.) As mentioned above, the principal 
challenge for numerical simulations is in approximating the fields p(x) and c(x) as they become 
unbounded. Our idea is to use a particle method for the evolution of p(x), so that its singularities 
manifest themselves as harmless particle aggregates. For the propagation of c(x), we use a second 
order implicit finite-difference scheme because of its simplicity and excellent stability properties. 
Even though c(x) also becomes unbounded analytically (it remains large but bounded for a given 
discretization), its weak logarithmic singularities do not cause any problems for sufficiently stable 
schemes. In general, one may use a different solver for propagation of the concentration field 
(including a particle-based solver as well). 

The particle-based description utilizes that the PDE ( [Ta| is a Fokker-Planck equation for the evolu- 
tion of the probability density of the bacteria. Suppose first, that we are given some time-dependent 
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Fi gure 1; A snapshot (t — 10) of the con- 
centration field c(x) illustrating multiple 
blow-ups in the parabolic model. The 
parameters are as follows: a = k = 1, 
x = 0.1, ]i = 0.005, M = 25. The spatial 
domain, 3.2 x 3.2, is discretized using the 
mesh size Ax = 0.05; the time step for 
o propagation of the concentration field is 

At = 0.1; 4096 particles are used to ap- 
proximate p(x) (see Section [5] for details 
of the numerical scheme). Some random 
initial data is used at t = 0. 

concentration field c(x, t) and want to solve equation |la| . Let xj"" 1 be i. i. d. random variables 
distributed according to the initial condition p(x, 0) at time t = and satisfying the following 
stochastic differential equations (with reflecting boundary at d Cl) for t > 0: 

dx| n) = X Vc(X t (,l) ,f) dt + ^/2fi dW t (n) . (3) 

Then, using the strong law of large numbers and that (la) is a Fokker-Planck equation for the SDEs 
l|3j, we can see that as N — >■ 00, the empirical probability densities, 

KM = ^Lt(x- x t n) )' ( 4 ) 

converge almost surely to p(x, t). (This convergence is in the sense of measures, point-wise in f; a 
stronger convergence can be established with more effort.) Thus we can simulate equations l|3| for 
the particles, and approximate the density p(x, t) using Pn{x, t) whenever required. 

Rigorously justifying this approach to the full K-S system (JTJ is more challenging because in 
the full system the concentration field c(x,t) is a functional of p(x,s), s < t, rather than some 
a priori prescribed field. This implies that equations (|3j become coupled and we end up with a 
system of stochastic interacting particles (with memory in the parabolic case). Some progress in 
this direction, nevertheless, has been made. Stevens derived the K-S equations as limit dynamics 
for interacting stochastic particle systems via smoothing and rescaling the interaction potentials in 
a particular fashion Il32l . Haskovec and Schmeiser proved convergence of their particle method 
under reasonable assumptions via analysis of the associated BBGKY hierarchy [22]. 




2. Numerical method 



As mentioned in the introduction, the principal idea behind our numerical scheme is to employ 
a particle method for the evolution of the particle density field p(x,t). The current implementa- 



tion utilizes an implicit second order finite-difference scheme for simulating equation I lb I and an 
explicit Euler-Maruyama scheme for the stochastic particle dynamics j3|. We discretize the compu- 
tational (rectangular) domain with grid size Ax and propagate the concentration field C« using the 
time step At. The particle density field P,j is reconstructed from the particle locations. We evolve 
the particles using adaptive time steps which may be smaller than At; this is needed for stability 
reasons (see below). 

The particle dynamics is simulated using the forward Euler-Maruyama scheme, 

xW(« + At) = xW(t) + *Vc(xW(*))At + ^/2^AtN(0,1), (5) 



3 




Figure 2: Adaptive time stepping for particle dynamics. In this ex- 
ample a blow-up of the concentration field occurs near a grid point 
with index ]. While some particular time step is appropriate for prop- 
agating the particles away from the blow-up (i), it may lead to over- 
shooting in its vicinity (ii, iii). Thus the particles incorrectly accumu- 
late near the grid points j ± 1 rather than One can avoid this phe- 
nomenon reducing the time step for the particle evolution whenever 
the gradient of the concentration field becomes too large. Practically, 
it is sufficient to adjust the time step so that the expected length of 
the particle jump does not exceed the mesh size Ax. 



where N(0, 1) is a standard Gaussian random variable with mean and variance 1. The gradient 
field Vc(x) is approximated in two steps. First, we construct the gradient fields CX and CY using 
the second order approximation, 
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Ci+y(t) - Q_y(t) , CYij = — Cy +1 (t) - Cy_!(f) . (6) 



Then we approximate Vc(ar) via bilinear interpolation using the values of CX and CY at the four 
nearest grid points. Because | Vc(x)| becomes unbounded (very large numerically) in the vicinity 
of blow-ups, it is essential to choose the time step At adaptively, see Figure [2] and its description 
for details. Consequently, if needed, we subdivide the timestep Af into smaller intervals Ar^ 
so that | Vc(xW) |AtW < Ax. (Each particle is simulated independently of others, i.e, At^"' are 
chosen for each particle individually.) A failure to satisfy this condition leads to overshooting in the 
vicinity of blow-ups, which in turn causes artificial oscillations and damping. Finally, we enforce 
the no-flux Neumann boundary conditions for p(x) by reflecting the particles escaping from the 
spatial domain back into it. 

The particle density field P;y(f) is reconstructed from the particle locations as explained in the de- 
scription to Figure [3] each particle contributes fractions of its weight to the four nearest grid points 
according to the bilinear interpolation rules. This approximation preserves the first moment of the 
corresponding distribution. It is important to use this kind of moment-preserving approximations, 
e.g., a simple bin-counting (assigning all weight of a given particle to the nearest grid point) is un- 
satisfactory. The reason is that this creates an artificial flux towards the grid points: the particles feel 
their own potentials which become artificially aligned to the nearest grid points. This flux is suf- 
ficient to pin the singularities to grid points and disrupt such phenomena as logarithmically-weak 
interaction of singularities with each other and with the boundary of the domain. 

The concentration field is propagated according to the implicit second order finite-difference 
scheme, 



a 

At 



Cij(t + At) - Q ; -(f)] = ^ D f C (* + At ) ~ + A + p <;(*)/ (7) 



where D^ 2 ' is the standard second difference operator. It is essential to use an implicit scheme for 
stability reasons because the concentration field acquires logarithmic singularities. This scheme 
is quite efficient: even though it requires solving a system of linear equations, the matrix of the 
finite difference operator is symmetric and banded, and may be Cholesky-factorized before the 
actual computations. The extra computational cost of solving this linear system is then alleviated 
by relaxation of the At ~ Ax 2 stability constraint of explicit schemes. 
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Figure 3: Schematic representation of the numerical scheme. 
The concentration field c(x) is sampled on a uniform grid (■), 
while the particles (-k) move without any restraints. The parti- 
cle density p(x) is reconstructed via bilinear interpolation: each 
particle contributes fractions of its mass to the four nearest grid 
points proportionally to the relative distances from the latter, 
e.g., a particle located at (x,y) within a unit-square grid cell 
contributes weights proportional to xy, x(l — y), y(l — x), and 
(1 — x)(l — y) to the respective grid points. 



The presented numerical method may be generalized to other similar systems of the form 

d t p(x,t) = V-(}i(x,p,c)Vp - pv{x,p,c)), (8a) 
a d t c(x,t) = V-(G(x)Vc) + F{x,p,c), (8b) 

and to more complicated domains. This makes it a useful tool for studying various aggregation 
phenomena for which the conventional finite-difference or finite-elements schemes are inefficient. 
From the numerical analysis point of view, we do not currently have any estimates of the conver- 
gence rates for this scheme. Partially, this problem is complicated due to the singular nature of 
solutions to the K-S equations. In particular, once the blow-ups are formed, one must interpret the 
K-S PDEs in some proper sense, and while our numerical scheme is able to propagate the solutions 
for all times (i.e., it regularizes the PDEs in some fashion), relation of such numerical solutions to 
a particular analytical regularization must be carefully investigated. On the other hand, we could 
speculate that our method exhibits the typical for this kind of schemes errors of the order 0(\/~Ki), 
0(1/VN), and 0(Ax 2 ) for regular solutions. Such analysis, however, is not the purpose of this 
paper and we verify our numerical results by "standard means," i.e., compare the results with 
simulations employing half the mesh size, double the number of particles, etc. 



3. Formation and interaction of singularities 

We now concentrate on the elliptic (a = 0) case, though some of the reasoning is equally applicable 
to the parabolic case as well. Generally, one can distinguish the soft and hard blow-ups: in the first 
case the fields p(x) and c(x) become unbounded but p{x) does not acquire atomic components, 
while in the second case the particle density field becomes a bona fide distribution. Here we con- 
centrate on the hard blow-ups which have a more apparent physical meaning: atomic components 
of the particle density correspond to accumulation of particles at some point locations on the given 
physical length-scale. 

It is convenient to study the elliptic model in a slightly different form, as a McKean-Vlasov 



system. First, we solve equation I lb I for c{x): 



c( x ) = - V d (x,y)p(y)dy = -V d *p(x). (9) 

Here V d {x,y) is the Green's function for A — k 2 in n, which is assumed to be d-dimensional. The 
entire-space Green's functions (the fundamental solutions) whose singular part is identical to those 
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in bounded domains are known explicitly: 

e -k\x-y\ 

Vi(x,y) - 
V 2 (x,y) -■ 
V 3 (x,y) -■ 



2k 
Ko(k\x 



In 

Q -k\x-y\ 

Att\x — y\ 



(d 
(d 
(d 



1) ; 

2) ; 

3) . 



(10a) 
(10b) 
(10c) 



Here Ko(-) is the modified Bessel function of second kind. Substituting expression |9) for the con- 
centration field c(x) into equation I la we get a closed integro-differential equation for the particle 
density field, 

d t p(x,t) = V-(jiVp + xpWd*p) = v- 



Sp 



Here the energy functional is given by 



£{p) = ]i 



f p(x)\np(x)dx + | JJ p(x)V i (x,y)p(y)dxdy. 
One of the features of the dynamics |Tl| is that the functional £{p) is non-increasing: 

M 



d£ 
dt 



n 



Sp 



p(x,t)dx < 0. 



(11) 



(12) 



(13) 



When p(x) becomes singular with respect to Lebesgue measure, the first (entropic) term in < |12[ | 
tends to infinity, which must be compensated by the second (interaction) term. If d = 1, the Green's 
function Vj (x) is bounded below and so is the interaction term, thus the blow-ups are not possible. 
This is no longer the case if d > 2 and so the blow-ups are permissible. A more delicate analysis is 
needed to understand how exactly they are formed. 

From standpoint of the stochastic process underlying the Fokker-Planck PDE (Ta| , the atomic 
components of p(x) correspond to particles aggregating at point locations. Thus the traps which do 
not allow the particles to escape must be created. Let us first investigate when a particle aggregate 
creates a trap for a single particle diffusing in its field. Suppose p(x) = MS(x), the corresponding 
concentration field is then proportional to the fundamental solution l |10) . Equation for the radial 
component of the location of a single particle diffusing in the field created by p is equivalent in law 
to the following SDE: 



dff 



X c'(r t ) + (d-1) 



dt 



2}i dW f . 



For d = 1 we get 



drt 



exp(-kr t )dt + ^2\i dW t . 



(14) 



(15) 



This equation is sufficiently regular and does not allow for existence of a trap for any value of M, so 
consistently with energy considerations, no blow-ups are possible in this case. For d > 3, a similar 
reasoning shows that a trap of an arbitrarily small mass may exist. So the interesting case is d = 2. 

For d = 2, using that in the leading order, as r — > 0, Ko(r) ~ — ln(r), we obtain the following 
equation: 

XM \ dt 



dr t 



1)1 dW f . 



(16) 



2n ) Tt 

This is a well-known Bessel process [27J, it's behavior near the boundary at r — depends strongly 
on the value of M. In particular, when 



M > 



Atzu 
X 



M* 



(17) 
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the origin is the so-called exit boundary, a trap, i.e., the particles reach it in finite time, and may never 
escape back to the domain r > 0. For our problem the implication is that the total mass of at least 
M* is needed for existence of a stable singularity in the particle density field. For M € (0, M*), 
the SDE ( [16) has the so-called regular boundary at r = 0. In this case the particles may reach the 
origin, but are not necessarily trapped there, i.e., may also leave it according to some rules which 
must be prescribed in addition to the SDE itself. This implies that the original K-S PDEs alone are 
not sufficient to describe the blow-up dynamics and the exchange of mass between the regular and 
singular components of the particle density when the mass of the singularity is less than M* . 

The critical mass M* is smallest mass such that a singularity does not shed its mass but only 
absorbs particles from the smooth component. It is, however, twice less than the mass required to 
create a singularity from smooth initial data. Indeed, consider the N-particle SDEs approximating 
the elliptic K-S equations: 

d ^" ) = -^^) ^ V d (xWxU)dt + ^dwW. (18) 

ly dX ' i=l, i^n 

To study the aggregation dynamics of the particle system, consider the variance of the empirical 
probability density of the bacteria Q, 

Rt will be referred to as the radius of the system. We concentrate on the problem in the entire space 
with k = 0; this corresponds to V%{x,y) = l/(2n) ln\x — y\. Analysis for the other values of k, or 
in the bounded domains is similar (although more technical) because the singularities of Vjj {x,y) 
are the same in the leading order. In this case, however, there exists an explicit closed equation 
governing the evolution of Rt'. 



The driving stochastic process Wf obeys the following SDE 




(20) 



1 N 

dWt = N^R t (^'^ " XfW ) ' dW ' 0, (21) 

and can be shown to be a Wiener process using the Levy characterization of Brownian motion. We 
can immediately see that in the limit as N — >■ oo, if 

M 7 L 8 Ht =Mc/ (22 ) 
X 

the evolution of the radius obeys a well-known deterministic equation, 

Rt = -^r, 7 = MM/M C -1). (23) 

Hence R| = Rq — jt , and for supercritical masses (M > M c ) the blow-up occurs at time T = Rq/J- 
It is independent of the initial distribution of mass except through Rq. If M = M c , after the time 
rescaling, equation |20| may be rewritten as 

dR t , = I| + dW f „ t' = f. (24) 
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This equation describes a critical Bessel process with entrance boundary at R = 0, in this case R t r a.s. 
never hits zero, although it visits its arbitrarily small neighborhood. In this limit the particle system 
remains stochastic even asN->oo and does not approximate the deterministic K-S equations. As 
t' — > oo, the particle system approaches blow-up infinitely often, but never actually coalesces into 
a single particle, cf (71 . 




Figure 4: Snapshots (f = 10000) of the concentration field c{x) — M/ (kL) 2 demonstrating sub- and supercriti- 
cal behaviors in the elliptic (a = 0) model. (The constant is subtracted because c(x) itself becomes unbounded 
in a finite domain as k — > 0.) The other parameters are as follows: k = 0.01, x = 0-1, ]i = 0.005; At = 0.01, 
Ax — 0.05, 4 • 10 3 particles. Initially, all particles are distributed randomly with a slight bias towards (0, 0). In 
the supercritical case (M = 0.35, right plot), a trap is created, while in the subcritical case (M = 0.34, left plot), 
all particles eventually spread uniformly over the entire domain. The theoretical critical mass for the elliptic 
model with k = in the entire space is M c = 2n/5 « 1.26. Notice that only a quarter of the critical mass is 
needed for a trap in the corner. 

Formation of singularities. Our first numerical experiment is designed to test how well our nu- 
merical method predicts the value of the critical mass M c required to develop a singularity from 
smooth initial data. Notice, first of all, that due to reflection principle for the problem with Neu- 
mann boundary conditions, singularity in a corner of the domain, O = (0, L) 2 , only requires a 
quarter of the critical mass. Since the singularity inside of the domain eventually migrates into 
a corner, we perform our experiment in the corner to begin with. We scatter all particles over fi 
with a slight bias towards (0,0) and vary their total mass observing whether the particles remain 
aggregated, or spread uniformly over the entire domain. The results are presented in Figure |4] 
For the specified values of the parameters, the calculated numerical value of the critical mass lies 
between 0.34 and 0.35, while the theoretical prediction for elliptic model with k = in the entire 
plane is M c /4 = 7r/10 — slightly smaller. The mismatch is a consequence of the simulation in a 
finite domain, the limit L — >■ oo, k — > is non-trivial because equation |Tb| with Neumann bound- 
ary conditions is not well-posed in Ci when a. = k = 0. In particular, the boundary of the domain 
far from the corner where the singularity is formed pulls the particles away from the singularity, 
effectively increasing the critical mass. 

Interesting phenomena occur when a singularity with mass M G (0, M* ) already exists in the 
system. In this case the trap absorbs the particles, though the latter may still escape back into the 
regular component of the particle density. The diffusion process underlying the K-S equations is 
not uniquely defined by its generator and additional rules for behavior of the particles at singular 
points must be specified. These rules are not inherently encoded in the K-S equations and are 
related to non-uniqueness of K-S regularizations. From the modeling perspective, the exchange of 
mass between the regular and singular components of the particle density in this regime strongly 
depends on particular details of the numerical method. Note, however, that if the initial particle 
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density is regular, the smallest possible mass for singularity is greater than M c = 2M* C , so this 
scenario does not occur. 




Figure 5: Snapshots of the concentration field c(x) illustrating the motion and interaction of singularities 
in the elliptic (a = 0) model. The other parameters are as follows: k = 1, x = 0.1, ]i = 0.005, M = 4; 
Af = 0.01, 4 • 10 3 particles. Initially, the particles are placed around four distinct locations; as time goes on, the 
singularities are formed, then they merge and travel towards the boundary of the domain, and finally stabilize 
in the corners. (Notice the scale change for the z-axis.) 



Interaction of singularities is illustrated in Figure [5] The singularities attract each other and are 
also attracted by the boundary of the domain. Equations governing their dynamics may be derived 
following [37. 38, 16J. If the density p(x) is purely atomic, we obtain 

Xi (t) = -X^-J2m j V 2 (xi,x j ). (25) 

Here m,-s are the masses of singularities located at Xj-s, ~Vi(x, y) is the Green's function. Note that 
unlike in the method of Haskovec and Schmeiser |21 1, in our method, the dynamics of singularities 
is not imposed by the numerical method explicitly, i.e., it is a natural consequence of the stochastic 
particle dynamics. 

Equations |25| are identical to equations <j3j of the K-S model with particle diffusivity y. set to 
zero. The naive reasoning is that once the traps are formed, the particles' diffusion is dominated 
by the (infinitely strong) drift and ceases to contribute into dynamics: the stochastic dynamics self- 
averages and the martingale component does not contribute into the mean drift. In order to verify 
how well our method approximates equations p5| , we perform a numerical experiment illustrated 
in Figure [6] In this experiment we create two singularities far from the boundaries in a sufficiently 
large domain, so that the interaction potential is well-approximated by the fundamental solution 
( jlObp . We run the K-S solver and track locations of these singularities. Finally, we compare them 
with locations of the particles evolving according to equation |25| from the same initial data. A 
perfect match indicates that the scheme is very successful in dealing with this kind of phenomena. 
Similarly, in a bounded domain, the particle interaction potential, V% {x\ — Xj), should be replaced by 
the Green's function which also encodes interaction with the boundary. In particular, this explains 
attraction of the singularities to the boundary and the corners of the domain. 
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3.1 



Figure 6: Interaction of singularities in the elliptic (ct = 0) 
model. The particle density p{x) is initialized with two 
delta-functions with masses 6.25 and 18.75. These singular- 
ities attract each other and eventually merge; their locations 
are plotted as functions of time t. The thin black line corre- 
sponds to the simulation of equation |25j, while the thick 
grey line — to simulation of the elliptic Keller-Segel model 
(the lines are barely distinguishable); At = 0.01, 4 ■ 10 5 par- 
ticles are used. The computational domain is chosen suffi- 
ciently large so that the boundary effects are negligible on 
the presented time scale. 
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3.2 



4. Discussion 

Our investigation of the K-S equations was motivated by the underlying stochastic particle dy- 
namics. A numerical method based on similar ideas has been recently presented by Haskovec 
and Schmeiser [21]. The greatest advantage of the H-S method is that the singularities are repre- 
sented as deterministic particles rather than as a tight cluster of "elemetary" stochastic particles. 
This, however, requires an explicit knowledge of evolution equations governing the interaction of 
singularities (at the very least the knowledge of Green's function in a given domain) and is only 
applicable to the elliptic K-S model. Moreover, the H-S method requires a direct simulation of an 
ensemble of interacting particles and thus its computational cost scales as N 2 with respect to the 
total number of particles. By combining the particle dynamics and the PDE dynamics of the con- 
centration field which mediates the particle interaction (the PIC ideas) we are able to significantly 
reduce this computational cost and deal with a much larger number of interacting particles. We 
can also treat both parabolic and elliptic models with equal ease. The ultimate numerical method 
should combine ideas presented in this work, the deterministic particles of the H-S method, and 
also high order conventional schemes in the regions where the solution remains regular. Such a 
multi-model implementation remains a subject for the future work. 

Even though the non-uniqueness issue is well-known, the details of how the exit-entrance con- 
ditions for the underlying diffusion process are related to the exchange of mass between the sin- 
gularity and the regular part are not yet well understood. From the modeling perspective, the 
implication is that a particular numerical scheme provides a regularization of some sort which 
affects the critical mass M* 

Another interesting class of stochastic particle models related to our studies of the K-S equations 
arises in the context of the so-called self -gravitating Brownian particles, see e.g., 1TTT1 11211311 . We 
suggest a model which fits into this class and bears an intrinsic relation to the K-S model. Consider 
an ensemble of particles in d dimensions characterized by their masses m n and evolving according 
to the stochastic differential equations 



Here No is the initial number of particles and M = J^m„ is their total mass. If initially m n = M/Nq, 
equation pet) is precisely < |18) . The interaction allows for particle collisions, thus equation pet) is 
only valid until the first collision, at which point the colliding particles coalesce into a single particle 
which acquires their combined mass. Dynamics is then restarted with the remaining particles. In 
the limit as No — > oo, whenever a particle accumulates anO(l) mass, its diff usivity becomes zero, 
i.e., it obeys the deterministic equation |25| which governs evolution of a point singularity in the 
K-S model. Therefore we conjecture that in a proper hydrodynamic limit this system is equivalent 
to the elliptic K-S model. Investigation of this system will be conducted elsewhere. 
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